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Abstract 

This  paper  compares  three  different  approaches  for  describing  the  cathode  catalyst  layer  of  a  PEMFC,  using  a  three-dimensional  CFD  model. 
The  three  catalyst  treatments  include:  a  thin-film  model,  a  discrete-catalyst  volume  model  and  an  agglomerate  model.  It  is  shown  that,  within  a 
single-phase  approach  using  physically  meaningful  parameters,  the  commonly  employed  thin-film  or  discrete-volume  descriptions  of  the  catalyst 
layer  do  not  show  the  significant  mass  transport  limitations  which  occur  at  higher  current  densities;  while  this  region  is  observed  using  a  catalyst 
agglomerate  approach.  Further,  an  in-depth  analysis  of  the  current  density  distributions  indicates  that  for  a  given  electrode  overpotential  the  thin- 
film  model  significantly  over-predicts  the  current  density,  compared  to  the  discrete  and  agglomerate  approaches.  The  thin-film  model  also  greatly 
exaggerates  the  variation  in  current  density  both  along  and  across  the  channel.  Finally,  the  agglomerate  model  predicts  noticeable  mass  transport 
losses  even  at  very  low  current  densities  despite  the  use  of  high  stoichiometric  ratios  and  thin-electrolyte  films. 

©  2008  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Fuel  cells,  particularly  the  polymer  electrolyte  membrane 
type,  have  long  been  touted  as  an  environmentally  friendly 
solution  for  automotive,  backup,  and  residential  power  needs. 
However,  despite  the  “green”  image,  fuel  cells  have  yet  to  be 
implemented  on  a  full  commercial  scale — due  in  a  large  part  to 
cost,  reliability,  and  lack  of  the  necessary  fuelling  infrastructure. 
The  optimization  of  the  PEMFC  related  to  cost  and  performance 
has  been  a  long  and  on-going  activity,  both  in  academia  and 
industry.  As  part  of  this  on-going  optimization  activity,  computa¬ 
tional  models  are  increasingly  being  applied  in  the  design  of  the 
catalyst,  porous  transport  layer  (PTL,  also  known  as  the  GDL), 
and  flow-field  plate;  especially  within  the  product  optimization 
and  research  prototyping  stages.  In  particular,  computational 
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fluid  dynamics  (CFD)-based  models  are  finding  increased  use 
due  in  part  to  their  ability  to  provide  detailed  spatial  resolu¬ 
tion  of  process  variables  such  as  species  concentration,  pressure, 
temperature,  electrical  potentials,  and  current  density  within  the 
individual  components  of  the  unit  cell.  These  variables  are  highly 
and  non-linearly  coupled  to  each  other  in  the  governing  conser¬ 
vation  equations  and  accordingly,  each  field  variable  influences 
and  is  influenced  by  the  electrochemical  reactions  occurring  in 
the  catalyst  layer.  Clearly,  it  can  be  expected  that  the  results  pro¬ 
duced  by  CFD  models  will  be  affected  by  the  extent  to  which 
the  physico-electro-chemical  processes  in  the  catalyst  layer  are 
considered.  More  specifically,  the  results  of  such  models  will 
depend  greatly  on  the  extent  to  which  gradients  in  species  con¬ 
centration,  temperature,  and  electrical  potentials  are  accounted 
for.  The  three  most  common  catalyst  layer  models  are — (i)  the 
agglomerate  model  [1,2],  (ii)  the  discrete- volume  model  [3-6], 
and  (iii)  the  thin-film  model  [7-9]. 

Among  the  three  approaches,  the  agglomerate  model  is  con¬ 
sidered  the  most  theoretically  detailed  as  it  attempts  to  include 
effects  due  to  the  catalyst’s  physical  structure.  In  this  approach, 
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Nomenclature 

aagg  effective  agglomerate  surface  area  (m2  m-3) 
apt  theoretical  platinum  loading  (m2  m-3) 
tfscale  scaling  factor  for  non-agglomerate  models 
C  local  species  concentration  (molm-3) 

D  diffusivity  of  dissolved  oxygen  in  the  electrolyte 
(m2  s”1) 

Dij  binary  diffusion  coefficient  (m2  s-1) 

Ed  Activation  energy  for  diffusion  (J  mol- 1 ) 

£th  Nemst  voltage  (V) 

Et  reaction  effectiveness  factor 

H  Henry  ’  s  constant  (Pa  m3  mol- 1 ) 

i  current  density  (A  m-3) 

i0  local  exchange  current  density  (A  m-2) 
kc  reaction  rate  constant  (s- 1 ) 

mpt  platinum  loading  (kg  m-2) 

Mw  molecular  weight  (kg  mol- 1 ) 

NCO  nominal  cathode  overpotential  (V) 

ragg  radius  of  the  agglomerate  particles  (m) 

rpt  platinum  particle  diameter  (m) 

S  volumetric  source  term 

t  thickness  (m) 

u  x-component  of  the  velocity  vector  (ms-1) 

v  y-component  of  the  velocity  vector  (ms-1) 

V\r  atomic  diffusion  volume  (m3) 

w  ^-component  of  the  velocity  vector  (ms-1) 

y  mass  fraction 

z  number  of  electrons  consumed  per  mole  of  reac¬ 

tant  (mol(e-)  mol(reactant)-1) 

Greek  letters 

ac  cathodic  transfer  coefficient 

<5  thickness  of  the  electrolyte  film  covering  the 

agglomerates  (m) 

£l  effective  platinum  surface  ratio 

rjact  activation  overpotential  (V) 

/ a  dynamic  viscosity  (kg  m-1  s-1) 

p  density  (kg  m- 3 ) 

ae  electronic  conductivity  (S  m- 1 ) 

cp  porosity  and  potential  (V) 

0agg  volume  fraction  of  the  electrolyte  phase 
0rh  relative  humidity 

<2>l  Theile’s  modulus 

Subscripts  and  superscripts 
agg  agglomerate 

cat  catalyst  layer  region 

e  electronic  or  electrode 

i  species  i 

op  operation  point 

p  protonic 

PTL  porous  transport  layer  region 

rxn  reaction 

sat  saturation  point 

th  theoretical  value 


Ionomer 


Fig.  1 .  A  typical  catalyst  agglomerate. 

the  catalyst  layer  is  considered  to  be  comprised  of  numerous 
agglomerates,  each  of  which  is  comprised  of  a  cluster  of  carbon 
black  particles  with  platinum  catalyst  dispersed  on  its  surface. 
The  Pt/C  catalyst  particles  are  held  together  by  a  polymeric 
electrolyte  material,  such  as  Nafion™.  The  catalyst  agglom¬ 
erates  are  generally  considered  to  be  of  one  of  the  following 
shapes — spherical,  cylindrical  or  slab;  with  the  pore  space  of 
the  catalyst  assumed  to  exist  among  the  percolating  network  of 
agglomerates  thereby  allowing  for  fluid  transport.  The  physical 
transport  processes  described  by  the  agglomerate  model  include 
the  gas-phase  transport  in  the  pore  space  of  the  catalyst  layer, 
dissolution  of  chemical  reactant  in  the  electrolyte  phase,  simul¬ 
taneous  diffusion  and  reaction  of  the  dissolved  reactant  within 
the  agglomerate,  ion  transport  in  the  electrolyte  phase,  and  elec¬ 
tron  conduction  via  the  carbon  black  particles.  A  representation 
of  a  single  spherical  agglomerate  can  seen  in  Fig.  1. 

The  physical  description  of  the  catalyst  layer  in  the  discrete- 
volume  model  is  very  similar  to  that  of  the  agglomerate  model 
with  a  few  exceptions.  The  discrete-volume  model  does  not  con¬ 
sider  the  formation  of  agglomerates  or  any  structure  beyond  the 
application  of  effective  medium  theory.  As  a  result  there  is  no 
consideration  of  the  dissolution  of  reactants  at  the  gas/Nafion™ 
interface,  nor  the  simultaneous  diffusion  and  reaction  within  the 
agglomerate  structure  itself.  In  one  published  variation  of  the 
discrete-volume  model  [10],  the  pore  space  of  the  catalyst  is 
considered  to  be  completely  filled  with  water. 

The  thin-film  model  is  very  different  from  the  discrete  or 
agglomerate  models;  in  this  approach,  the  catalyst  layer  is  not 
explicity  described  but  is  treated  as  a  boundary  condition  at  inter¬ 
face  of  the  PTL  and  the  membrane.  The  treatment  of  the  catalyst 
layer  as  an  interface  means  that  it  does  not  account  for  any  trans¬ 
port  or  resistance  in  the  physical  structure  of  the  catalyst  layer. 

In  this  paper,  a  half-cell  PEMFC  cathode  model,  developed 
in  our  research  group,  is  applied  to  compare  the  predictions 
made  using  the  three  different  catalyst  layer  models.  The  basis 
of  comparison  among  the  different  models  is  made  in  terms  of 
their  predicted  polarizations  and  spatial  distribution  of  current 
densities.  The  results  are  interpreted  in  the  context  of  two  key 
factors  affecting  the  current  density  and  its  distribution — oxygen 
concentration  and  overpotential. 

2.  Computational  domain 

The  fuel  cell  considered  in  this  paper  represents  a  single  chan¬ 
nel  of  a  parallel  flow  field,  and  as  such,  the  model  domain  is 
greatly  simplified  compared  to  a  serpentine  flow-field  in  which 
convection  in  the  PTL  [11]  must  be  considered  due  to  the 
in-plane  pressure  gradients  [12-14],  A  typical  computational 
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Fig.  2.  Three-dimensional  domain  of  a  cathode  with  straight  channels. 


Table  1 

Fuel  cell  geometry 


Parameter 

Value 

Unit 

Channel  length,  L 

80 

mm 

Channel  width,  ttiCh 

1 

mm 

Channel  height, 

1 

mm 

Channel  hydraulic  diameter,  D h 

1 

mm 

Land  width,  ini 

1 

mm 

PTL  thickness,  ipri. 

300 

(urn 

PTL  porosity,  0ptl 

0.4 

Catalyst  thickness,  tc at 

12.9 

jjum 

Catalyst  porosity,  0cat 

0.1 

exception  of  the  catalyst  layer  thickness  for  the  thin-film  model 
which  is  zero — consistent  with  the  assumptions  of  the  approach. 

3.  Transport  phenomena 

The  cathode  model  is  single  phase,  non-isothermal  and 
accounts  for  electron/reactant  transport  in  the  PTL  and  pro¬ 
ton/electron/reactant  transport  in  the  catalyst  layer.  The  transport 
equations  are  solved  using  the  commercial  CFD  solver  Fluent 
6.1  with  the  different  catalyst  models  implemented  through  the 
use  of  custom- written  user  defined  functions  that  effectively  add 
source  and  sink  terms  within  the  catalyst  layer,  enforce  specific 
boundary  conditions,  and  implement  specific  material  proper¬ 
ties.  Although  the  model  is  non-isothermal  in  that  it  accounts 
for  heat  generation  due  to,  both,  reaction  and  activations  losses 
it  was  run  isothermally  for  this  investigation.  It  is  important  to 
emphasize  that  the  accuracy  of  thermal  boundary  conditions, 
in  particular  the  plate/coolant  channels,  and  the  material  prop¬ 
erties  of  both  the  MEA  and  plate  are  critical  in  achieving  the 
correct  temperature  distribution  through  the  MEA  and  are,  for 
the  moment,  outside  the  scope  of  this  work.  All  transport  proper¬ 
ties  within  the  porous  media  are  assumed  to  be  isotropic,  and  the 
Bruggeman  relation  is  used  to  estimate  the  effective  mass  dif- 
fusivities  in  the  porous  electrodes.  Further  details  of  the  model 
are  given  in  Ref.  [13], 

3.1.  Gas  flow  channels 

The  flow  within  the  gas  supply  channels  was  modelled  as  an 
ideal  gas  flow  with  the  following  governing  equations. 


domain  used  in  this  work  is  shown  in  Fig.  2;  the  length  of  the 
cell  used  in  this  study  was  80  mm  with  a  small  inlet  and  out¬ 
let  added  upstream  and  downstream  of  the  active  area  to  port 
the  gases  to  the  flow  channel.  The  domain  represents  the  cath¬ 
ode  flow-field,  a  PTL,  and  a  catalyst  layer.  The  computational 
domain  was  reduced  to  only  one-half  a  landing  and  one-half 
a  channel,  as  is  shown  in  the  cross-section  of  Fig  3,  since  the 
geometry  and  boundary  conditions  are  symmetric  about  the  line 
x / 1\  =  0  (this  is  reasonable  if  convective  effects  -  which  can 
occur  in  straight  channels  due  to  poor  plate  manifolding  -  are 
neglected).  Details  of  the  geometry  are  given  in  Table  1 ,  with  the 


Fig.  3.  Cross-section  of  the  computational  domain. 


3.1.1.  Conservation  of  mass 

^-{pu)  +  ^-(pv)  +  ^-{pw)  =  0  (1) 

ox  3 y  dz 

where  p  is  the  mixture  density. 

The  density  of  the  fluid  mixture  is  determined  using  an  ideal- 
gas  mixing  law: 

Pmixture  =  -  (2) 

RTYpi/Mw,d 

where  P0 p  is  the  operating  pressure,  MW)1-  is  the  molecular  weight 
of  species  i,  T  is  the  temperature,  R  is  the  universal  gas  constant, 
and  yi  is  the  mass  fraction  of  species  i. 


3.1.2.  Conservation  of  momentum 


3  (pu) 

d(pu) 

d(pu) 

dp  2  ,  c 

(3) 

dx  dy  11 

dz  - 

-  ~z~  +  p  V  u  +  Smom* 
dx 

3  Ipv) 

d(pv) 

d(pv)  _ 

dp  o 

— - b  V  V  +  Smorm, 

dy 

(4) 

U  «  r  V  ~  r  U) 

ox  dy 

dz 

d(pw)  . 

d(pw)  . 

d(pw) 

dp  7 

,z  (5) 

dx  dy 

v  dz 

~  g-  +  tt  V  W  +  Smon 

where  5m0mx ,  Smoniy ,  and  Smom-  are  momentum  source  terms  for 
the  x,  y,  z  momentum  equations,  respectively,  p  is  the  dynamic 
viscosity  of  the  fluid  mixture  and  is  determined  using  a  semi- 
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empirical  formula  [15]: 

»*= E w 


In  the  present  model,  the  superficial  velocity  was  used;  how¬ 
ever,  the  mass  flow  rate  entering  the  porous  media  is  determined 
using  the  superficial  velocity  and  as  such  the  total  pressure 
drop  across  the  media  was  the  same  for  both  formulations 
[16]. 


with  <pij  determined  by 

[i  4  2 

0ij  J1/2 

and  the  mole  fraction,  x*,  as 
yi/Mwj 

Ypi/M^j) 

i= 1 


(7) 


(8) 


3.1.3.  Species  transport 

For  the  purposes  of  this  model,  the  dilute  approximation 
for  species  transport  was  not  considered;  rather,  a  full  multi- 
component  treatment  of  species  diffusion  was  used.  The  general 
form  of  the  equation  describing  the  species  transport  is  [16]: 

V  •  (pvy,)  =  -V  •  j,  (9) 

where  J,  is  the  diffusive  mass  flux  vector  for  species  i  and  was 
described  by  the  Maxwell-Stefan  equation  for  multicomponent 
mixtures: 

N- 1 

J,  •  -5>DyV>V  (10) 


The  Maxwell-Stefan  diffusion  coefficients  can  be  approx¬ 
imated  with  the  binary  coefficients  [17]  and  were  determined 
using  the  Fuller,  Shettler,  and  Giddings  empirical  correlation 
[18]: 


T^75((l /MWii)  +  (1  /MWJ))1/2 

P((AVM)1/3 -K£fcVkj)1/3)  * 


(11) 


where  V4,;  and  V4j  are  the  atomic  mass  volumes  for  species  i 
and  j,  respectively;  these  values  are  tabulated  in  Cussler  [18]  for 
several  common  compounds. 


3.2.2.  Conservation  of  momentum 

The  porous  media  in  this  study  was  modelled  through 
the  inclusion  of  a  momentum  source  term  within  the  general 
momentum  equations;  this  source  term  consisted  of  two  sub¬ 
terms — a  viscous  loss  and  an  inertial  loss: 

Smom;  =  —  PijHvj  +  y,  Cjj  ^  pvm  ag  (14) 

where  C  and  D  are  diagonal  matrices  that  contain  the  inertial 
resistance  factor  and  inverse  of  the  permeability,  respectively. 

A  Darcy  flow  regime  is  assumed  for  flow  in  the  porous  media, 
such  that  the  non-linear  term  is  negligible  and  Eq.  (14)  reduces 
to 


where  k  is  the  permeability  of  the  porous  media. 


3.2.3.  Species  transport 

The  species  transport  equation  in  the  porous  media  was  as 
follows,  with  a  modification  to  account  for  the  porosity  [19]: 

V  ■  {pvyd  =  -V  ■  (</.ptlJ,)  (16) 

where  v  is  the  superficial  velocity  (recalling  the  relationship  with 
the  physical  velocity  in  Eq.  (13)).  The  diffusive  flux  vector,  J,-, 
is  corrected  using  an  effective  binary  diffusivity  value  which 
reflects  the  existence  of  the  porous  media: 

N- 1 

J;  =  -ypD^Vyj  (17) 

j=l 

where  De g-  is  determined  by  a  Bruggemann-type  correction  [20] : 

Deft  =  Dijf1^ l  (18) 


3.2.  Porous  transport  layer 

The  PTL  is  a  porous  media  and  as  a  result,  the  governing 
equations  presented  for  the  flow  channel  still  apply;  but,  with 
modifications  that  account  for  the  porous  effects. 

3.2.1.  Conservation  of  mass 

V  •  (pUDarcy)  =  0  (12) 

where  roarcy  is  the  Darcy  or  superficial  velocity  and  is  based  on 
the  volumetric  flow  rate  of  the  porous  media  [16].  The  superficial 
velocity  is  related  to  the  physical  or  intrinsic  velocity  of  the 
porous  media  by  the  Dupuit-Forcheimer  relationship  [19]: 

I’Darcy  =  ^PTL^physical  (13) 

where  ^physical  is  the  local  volume  average  of  the  fluid  velocity. 


3.2.4.  Conservation  of  charge 

The  transport  of  electrons  through  the  solid  phase  of  the 
porous  medium  was  described  using  Ohm’s  law  [21]: 

-V  •  (ere  V 0e)  =  0  (19) 

where  ere  is  the  specific  electronic  conductivity  of  the  solid 
material  in  the  PTL  and  <fic  is  the  electronic  potential. 

3.3.  Catalyst  layer 

The  catalyst  layer  was  modelled  as  a  porous  media,  simi¬ 
lar  to  the  PTL;  therefore,  the  equations  governing  the  physics 
of  the  PTL  are  the  same  as  those  within  the  catalyst  with  the 
exception  that  the  catalyst  region  includes  an  extra  equation 
that  governs  the  transport  of  protons  and  an  appropriate  set  of 
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sources  and  sinks  that  account  for  the  processes  of  generation 
and  consumption  based  on  the  cathode  half-cell  reaction. 

3.3.1.  Conservation  of  mass 

v  ■  (PV)  =  Sm ass  (20) 

where  Smass  is  a  source  of  mass  per  unit  volume: 

Smass  =  Sh20  +  S02  (21) 

with  .S'h20  and  Sq2  as  the  source  terms  for  water  and  oxygen  per 
unit  volume,  respectively. 

3.3.2.  Conservation  of  momentum 

Recalling  the  Darcy-like  form  of  Eq.  (15),  the  momentum 
equations  for  the  catalyst  layer  are 

V  ■  ( pvv )  =  —V p  +  i u.V2u  —  +  Smom„.  (22) 

k 


concentration  gradient.  These  two  effects,  electro-osmotic  drag 
and  back-diffusion,  are  generally  combined  into  a  single  expres¬ 
sion  known  as  the  net  drag  coefficient  of  water  [24,25,27,29]. 
This  net  drag  coefficient  was  measured  by  Choi  et  al.  [30]  and 
was  then  correlated  with  current  density  and  cathode  potential 
by  Sun  et  al.  [1]: 

{1.0,  NCO<  0.25  V 

46  x  NCO2  -  31.52 

xNCO  +  5.7,  0.25  V  >  NCO  <  0.35  V 

0.3,  NCO  >  0.35  V 

(26) 

where  NCO  is  the  nominal  cathode  overpotential  and  is  outlined 
in  Eq.  (37). 

The  overall  contribution  of  the  water  flux  across  the  mem¬ 
brane  is  described  by 


where  .SVnomrxr,  is  the  source  term  for  the  ith  momentum  equation 
in  x,  y,  z  which  is  included  to  account  for  the  effect  of  the  added 
mass  resulting  from  the  electrochemical  reaction: 

5momrxn  =  SmassVi  (23) 

3.3.3.  Species  transport 

The  transport  of  species  is  governed  by  the  following,  with 
the  addition  of  a  source  term  to  account  for  the  consump¬ 
tion/generation  of  reactants/products  by  the  electrochemical 
reaction: 


2«Mh2o  . 

5H20flux  =  — -  i 

zh2oF 

The  total  source  of  water  is  then 
Sh20  =  ^H2Orx„  +  SH2Oaux 


(27) 


(28) 


3.3.4.  Conservation  of  charge 

Electron  transport.  Electrons  in  the  catalyst  are  transported 
through  the  solid  phase  (platinum  and  the  graphite  or  carbon 
supporting  matrix)  by  means  of  conduction: 


V  ■  ( pvyi )  =  -  V  •  (0catJ f)  +  Si 


(24)  -V  ■  (crecatV0e)  =  Se 


(29) 


where  5,  is  the  source  of  species  i  per  unit  volume. 

The  source  terms  for  each  species  were  determined  based 
on  the  cathode  half-cell  reaction  with  the  following  equation 
representing  the  amount  of  species  i  (oxygen  or  water)  which 
was  consumed/generated  by  the  electrochemical  reaction: 


where  Zi  is  the  number  of  electrons  consumed  per  mole  gener¬ 
ated/consumed  of  species  i,  F  is  Faraday’s  constant,  and  i  is  the 
local  source  of  current  per  unit  volume. 

In  addition  to  the  water  generated  by  the  electrochem¬ 
ical  reaction,  there  is  also  a  movement  of  water  between 
the  anode  and  the  cathode.  Typically,  this  interchange  of 
water  is  considered  to  be  the  result  of  three  separate 
mechanisms — electro-osmotic  drag  due  to  the  potential  gradient 
in  the  charge  carrier,  back-diffusion  resulting  from  a  concentra¬ 
tion  gradient  in  water,  and  convection  due  to  pressure  gradients 
that  arise  from  membrane  swelling  (capillary  pressure  and  elas¬ 
tic  stresses)  [22-24] .  The  transport  of  water  across  the  membrane 
due  to  the  mechanism  of  convection  is  not  considered  and  is 
referred  for  future  work.  Generally,  the  movement  of  water  due 
to  electro-osmotic  drag  is  described  through  the  use  of  a  drag 
coefficient  that  describes  the  number  of  water  molecules  trans¬ 
ported  per  charge  carrier  (H+),  this  value  ranges  from  0.1  to  2.5 
in  the  literature  [25-28],  The  back-diffusion  of  water  from  the 
cathode  to  the  anode  is  expressed  as  a  function  of  the  local  water 


where  Se  is  the  local  source  of  electric  current  per  unit  volume 
and  within  the  catalyst  layer  the  specific  electronic  conductiv¬ 
ity  is  a  function  of  the  catalyst  porosity  and  Nation™  volume 
fraction  in  the  agglomerate: 

°ecat  =  °e[(l  -  tfcatXl  -  0agg)]''5  (30) 

Proton  transport.  Protons  generated  in  the  anode  reaction  are 
transported  by  the  solid  polymer  membrane  from  the  anode  to 
the  cathode  and  are  transported  in  the  catalyst  layers  via  the 
Nafion™.  This  process  can  be  described  by 

-V  ■  (CTpV0p)  =  Sp  (31) 

where  5p  is  the  local  source  of  protonic  current  per  unit  volume, 
0p  is  the  protonic  phase  potential,  and  crp  is  the  specific  protonic 
conductivity. 

The  conductivity  of  the  Nafion™  phase  that  exists  within  the 
catalyst  is  dependent  on  the  local  water  content  and  is  assumed 
to  follow  the  relationship  reported  by  Springer  et  al.  [31]  (this 
relationship  will  vary  with  different  ionomers  types): 

Op=100(5.139e-3A.— 3.26e-3)exp  [l268  (32) 

where  A.  is  the  local  water  content  of  the  Nafion™  and  is  a 
function  of  the  local  relative  humidity.  This  equation  is  valid 
for  relative  humidities  less  than  or  equal  to  unity — above  unity 
Schroeder’s  Paradox  occurs  (given  the  single-phase  nature  of 
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this  work-  -Shroeder’s  Paradox  is  currently  neglected).  The  rela¬ 
tionship  between  local  humidity  and  membrane  water  content 
was  described  by  Hinatsu  et  al.  [32]: 

A-  =  0.3  +  1O.80rh  -  164h  +  14. 1 4H  (33) 


where  </>rh  is  the  relative  humidity  and  was  determined  by 


^t20,sa 


(34) 


with  Ph2o  as  the  local  partial  pressure  of  water  and  Pf^o.sat  as  the 
local  saturation  pressure  of  water — determined  by  a  regression 
of  the  data  in  Ref.  [33]: 


PH2o,sat  =  1.26837  x  10“3r4  -  1. 49827 T3 

+6.70916  x  102  r2  -  1.34832  x  105r 
+  1.02503  x  107  (35) 


3.4.  Activation  overpotential  and  cathode  potential 


The  activation  overpotential  is  a  measure  of  the  voltage 
losses  that  result  from  driving  the  electrochemical  reaction  at  the 
electrode  [34],  Through  consideration  of  the  phase  potentials, 
the  activation  overpotential  can  be  evaluated  via  the  following 
[1,35]: 

Oact  =  0e,local  —  0p, local  —  0ref  (36) 

where  0e, local  is  the  local  electronic  potential,  <£p, local  is  the  local 
protonic  potential,  and  0rer  is  the  desired  reference  electrode. 

With  respect  to  this  model,  the  reference  electrode  has 
been  chosen  as  the  oxygen  electrode  (similar  to  that  of  Sun 
et  al.  [1,36]);  as  such  the  reference  electrode  potential  is 
zero. 

The  cathode  potential  or  nominal  cathode  potential  (NCO) 
is  the  difference  between  the  electronic  phase  potential  at  the 
collector  plate  (bi-polar  plate)  and  the  protonic  phase  potential  at 
the  catalyst/membrane  interface.  This  term  is  a  description  of  the 
total  losses  in  the  cathode — activation  losses,  ohmic  losses  due 
to  electron  transport  in  the  catalyst  and  PTL,  ohmic  losses  due 
to  protonic  transport  in  the  catalyst,  and  mass-transport  losses. 
It  follows  that  the  operating  cathode  voltage  can  be  described 
by 

Kathode  =  Eth  -  NCO  (37) 


where  £+  is  the  Nemst  voltage  [34]: 


£th  =  E°  +  —  In 


(38) 


where  E°  is  the  open  circuit  voltage  at  the  specific  operating 
temperature  and  an2g ,  «o2g ,  «h2Oi  are  the  activities  of  gaseous 
hydrogen,  gaseous  oxygen,  and  liquid  water,  respectively. 

For  the  purposes  of  determining  the  open  circuit  voltage;  the 
total  pressure  on  the  anode  is  assumed  to  be  1  atm  and  have  a 
humidification  level  of  100%. 


3.5.  Approaches  to  catalyst  layer  modelling 

The  following  section  outlays  the  methodologies  and  equa¬ 
tions  for  each  catalyst  modelling  approach:  agglomerate, 
discrete-volume,  and  thin-film  interface. 


3.5.1.  Agglomerate  catalyst  model 

In  this  approach  the  transport  processes  within  the  catalyst 
are  assumed  to  be  described  by  the  following: 

•  Dissolution  of  oxygen  across  the  gas-Nafion™  interface. 

•  Diffusion  of  dissolved  oxygen  within  the  Nafion™  elec¬ 
trolyte. 

•  Diffusion  of  dissolved  oxygen  within  the  agglomerate  struc¬ 
ture. 


The  agglomerate  catalyst  model  is  described  by  the  standard 
Butler- Volmer  kinetics  [35]: 

.  eff.  C02  [  f  acF  \  /(I  -ac)F  \1 

[e*P  (f*r  -eXP 

(39) 


following  a  detailed  re-arrangement,  with  consideration  for  the 
mechanisms  of  oxygen  transport  within  the  catalyst,  the  govern¬ 
ing  kinetics  take  the  following  form  [1]: 


zFPQ2  (  1 _ (ragg  +  8)8 \  1 

H  V  £rkc(l  -  0cat)  flagg'+ggO  ) 


(40) 


where  Er  is  an  effectiveness  factor  of  the  spherical  agglomerate, 
z  is  the  number  of  electrons  involved  in  the  reaction  per  mole  of 
reactant,  H  is  Henry’s  Constant,  aagg  is  the  effective  agglomerate 
surface  area,  ragg  is  the  agglomerate  radius,  <5  is  the  thickness  of 
the  Nafion™  film,  D  is  the  diffusivity  of  the  dissolved  oxygen 
in  Nafion™,  and  kc  is  a  reaction  rate  constant. 

The  effectiveness  factor  of  the  spherical  agglomerate  is  also 
a  measure  of  the  effectiveness  of  the  electrode  reaction  and  can 
be  found  via: 


,-JA 


<Z>l  \tanh(3^L)  30l  / 


(41) 


with  0[  ,  Theile’s  modulus  for  chemical  reactions,  as 


(42) 


where  Dcfr  is  the  effective  diffusivity  of  dissolved  oxygen  in 
Nafion™  within  the  agglomerate  structures  and  is  determined 
using  a  Bruggemann-type  correction: 

OelT  =  D0lg5g  (43) 

with  0agg  as  the  volume  fraction  of  Nafion™  in  the  porous  media 
and  D  as  the  diffusivity  of  dissolved  oxygen  in  Nafion™  [37]: 


D  =  4.38  x  1(TS  exp  — 


(44) 
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Ed  is  the  activation  energy  for  diffusion  and  was  determined 
from  the  data  presented  by  Parthasarathy  et  al.  [37]. 

Henry’s  constant  is  determined  using  the  relationship  of  [38]: 


7T  =  7T  exP - 


(45) 


where  A  Gs  is  the  free  energy  of  dissolution  and  is  approximately 
-5209.69  J  mol- 1 . 

The  reaction  rate  constant,  kc,  is 


< 

ZF{  1  -  0cat)l 


Table  2 

Operating  conditions 


Parameter 

Value 

Unit 

Stoichiometric  ratio,  A 

5.0 

Inlet  relative  humidity,  0rh 

0.5 

Inlet  temperature,  T 

50 

°C 

Operating  pressure,  P 

1.5 

atm 

is  the  local  gas  concentration  of  oxygen;  this  is  a  major  differing 
point  with  the  agglomerate  model  which  instead  considers  the 
concentration  of  dissolved  oxygen. 


-  exp  ^ (1  RaTc)F  t?act)  ]  (46) 

where  ajlf  is  the  effective  platinum  surface  area  per  unit  catalyst 
layer  volume,  i0  is  the  exchange  current  density  at  temperature, 
T,  Cr(f2  is  the  reference  oxygen  concentration,  and  rjaci  is  the 
local  activation  overpotential. 

The  exchange  current  density,  i0,  was  calculated  using  the  ref¬ 
erence  exchange  current  density,  corrected  for  local  temperature 
(this  form  was  used  in  each  catalyst  model): 


i0  =  i0  exp - 


(47) 


3.5.3.  Thin-interface  catalyst  model 

The  thin-interface  approach  treats  the  catalyst  in  a  manner 
that  it  is  considered  as  only  a  layer  of  source  terms  for  reactants, 
energy,  and  electrons  [39].  This  is  achieved  by  considering  the 
catalyst  as  having  only  a  single  control  volume  in  the  thick¬ 
ness  direction  and  neglecting  the  heat  transfer,  proton  transport, 
electron  transport,  reactant  transport,  and  activation  overpoten¬ 
tial  distributions  that  occur  within  the  catalyst.  This  method  is 
governed  by  the  same  form  of  the  Butler- Volmer  kinetics  as 
shown  in  Eq.  (50). 

3.6.  Operating  conditions 


where  i™*  is  the  reference  exchange  current  density  at  the  refer¬ 
ence  temperature,  T0. 

aptf  is  a  function  of  the  specific  surface  area  per  unit  catalyst 
volume,  apt  and  an  effectiveness  factor,  £l,  which  accounts  for 
the  portions  of  the  catalyst  that  are  unable  to  meet  the  require¬ 
ments  for  electrochemical  reaction  (i.e.  electrically,  protonically, 
and/or  gas-phase  isolated): 

apf  =  £L«Pt  (48) 


apt  was  determined  from: 
3mpt 


(49) 


TptPPtlcat 

where  mpt,  ppt,  and  tcat  are  the  platinum  loading,  density  of 
platinum,  and  catalyst  layer  thickness,  respectively. 


3.5.2.  Discrete-catalyst  model 

Within  this  method,  the  catalyst  is  considered  to  be  governed 
by  the  standard  Butler- Volmer  type  kinetics,  as  in  the  agglom¬ 
erate  model,  however  this  method  gives  no  consideration  for  the 
transport  of  oxygen  through  the  electrolyte  phase.  The  model 
does  consider  the  resistance  posed  by  the  layer  for  reactant  and 
charged  species  transport. 

The  standard  Butler- Volmer  kinetics  are 

.  Co2  [  (  «c  F  \  /(I  -ac)F  Yj 

*— flscale*o  £,ref  exP  (  7act  I  exP  (  Oacl J 

(50) 

where  ascaie  is  a  scale  factor  used  to  adjust  for  the  active  area 
of  the  catalyst  in  order  to  calibrate  the  cathode  polarization 
curves  for  each  of  the  non-agglomerate  approaches  and  Cq2 


The  operating  conditions  employed  in  this  work  are  given 
in  Table  2  while  the  catalyst  parameters  are  given  in  Table  3. 
The  stoichiometric  ratio  was  deliberately  set  high  (A.  =  5.0)  in 
order  to  minimize  the  mass  transport  effects  which  occur  due 
to  reactant  consumption.  If  significant  differences  between  the 
approaches  are  noted  under  these  conditions,  they  are  certain 
to  be  of  even  more  importance  under  realistic  flow  rates.  Also 
of  note  is  the  fact  that  the  inlet  relative  humidity  was  held  at 
50%  in  order  to  prevent  the  formation  of  liquid  water  in  the  PTL 


Table  3 

Catalyst  and  PTL  properties 


Parameter 

Value 

Unit 

Reference 

Platinum  loading,  mpt 

4 

gm"2 

[13] 

Platinum  radius,  rpt 

1.5 

nm 

[1] 

Agglomerate  radius,  ragg 

1 

(jum 

[1] 

Eff.  agglomerate  area,  aa gg 

3.6  x  105 

m2  m-3 

[1] 

Ref.  Exch.  current  density,  i[f 

£celI>0.8V 

3.85  x  10“4 

A  cm-2 

[38] 

£cetl<0.8V 

1.5  x  10“2 

A  cm”2 

Activation  energy,  E^i 

-Been  >  0.8  V 

76.5  x  103 

Jmol  1 

[38] 

Ecc„  <  0.8  V 

27.7  x  103 

J  mol-1 

Charge  transfer,  ac 

Eceii  >  0.8  V 

Eceii  <  0.8  V 

1.0 

0.55 

[38] 

Ref  O2  concentration,  Co2,rcf 

0.85 

molm  3 

[38] 

Henry’s  constant,  H 

0.2685  x  105 

Pam3  mol-1 

[37] 

Effective  Pt  surface  ratio,  £l 

0.75 

[41] 

Electrolyte  fraction,  0agg 

0.5 

[42] 

PTL  conductivity,  <tptl 

100 

Sm"1 

[13] 
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and  channels,  as  the  physics  of  liquid  water  transport  are  not 
included  in  this  model.  An  inlet  value  of  </>rh  =  50%  was  found 
adequate  to  keep  the  relative  humidity  everywhere  below  100%, 
and  this  is  supported  experimentally  by  recent  neutron  imaging 
data  [40], 

3. 7.  Boundary  conditions 

A  fully  developed  laminar  flow  profile  is  assumed  at  the  chan¬ 
nel  inlet  with  the  properties  given  in  Table  2.  The  electronic 
phase  potential  is  set  to  zero  on  the  landing  surface  and  a  zero 
flux  of  electrons  is  imposed  on  the  membrane/catalyst  interface. 
The  ionic  phase  potential  is  set  to  the  desired  load  point  at  the  cat¬ 
alyst/membrane  interface  and  a  zero  flux  of  protons  is  enforced 
at  the  catalyst/PTL  boundary.  All  equations  employ  symmetry 
on  the  x/Dh  =  —  1  surface  and  the  x/D^  =  Osurface.  The  outlet 
condition  assumes  a  zero  gradient  in  the  flow  direction  for  all 
variables,  with  the  exception  of  pressure. 

4.  Results 


Fig.  5  shows  the  predicted  current  density  distribution  nor¬ 
malized  by  the  load  point  average  current  at  a  location  very 
near  the  inlet.  The  results  are  in  good  agreement  with  the  two- 
dimensional  predictions  of  Sun  et  al.  [  1  ]  and  aid  in  understanding 
the  limiting  factor  in  current  production  for  the  chosen  load 
points.  At  each  point,  the  local  reaction  rate  is  dependant  upon, 
both,  the  oxygen  concentration  at  the  reaction  site,  and  the  local 
overpotential;  which  is  a  function  of  the  local  ionic  end  electronic 
potentials. 

At  low  currents  (load  point  A)  there  is  relatively  little  oxy¬ 
gen  consumption,  and  as  such  the  reaction  rate  is  limited  by  the 
transport  of  electrons  through  the  PTL  which  results  in  the  max¬ 
imum  current  production  occurring  under  the  land  area  where 
electrons  are  most  accessible.  As  the  current  is  increased  (load 
point  B);  both,  electron  transport  and  oxygen  transport  become 
important  such  that  the  location  of  maximum  current  production 
shifts  to  the  interface  between  the  land  and  the  channel.  Finally, 
at  high  currents  (load  point  C),  the  transport  of  oxygen  limits  the 
reaction,  and  the  maximum  current  production  is  now  occurring 
under  the  channel  where  the  oxygen  is  most  abundant. 

A  significant  observation  here  is  that  the  effect  of  the  mass 
transport  limitations  are  clearly  evident  in  the  case  of  load  point 
B  even  though  this  load  point  is  firmly  located  in  the  linear,  or 
‘classic’  ohmic,  region  of  the  polarization  curve.  In  hindsight, 
this  should  not  be  an  unexpected  effect;  polarizations  curves  are 
concave  up  at  low  current  and  concave  down  at  high  current  and 
there  must  be  an  inflection  point  somewhere  in  between.  This 
inflection  point  is  the  location  where  mass  transport  losses  actu¬ 
ally  begin  to  become  significant  and  is  usually  quite  subtle  as  the 
dominating  loss  is  transitioning  from  one  mechanism  to  another. 

We  now  turn  our  attention  to  the  three-dimensional  effects 
that  occur  along  the  channel.  For  both  load  points  A  and  C 
(not  shown),  the  current  density  distributions  are  very  similar 
to  those  in  Fig.  5  with  a  decreasing  magnitude  as  the  distance 
along  the  channel  increases.  However,  this  is  not  the  case  for 
load  point  B;  Fig.  6  depicts  the  current  density  distribution  for 
load  point  B  at  three  locations  along  the  channel.  Near  the  inlet, 
the  maximum  current  density  is  located  near  the  land/channel 


4.1.  Catalyst  agglomerate  model 

The  catalyst  agglomerate  model  is  the  most  detailed  of  the 
catalyst  models  described  in  this  paper  in  that  it  accounts  for 
more  physical  processes  and  a  better  description  of  the  actual 
morphology  of  the  catalyst  layer.  As  such,  it  is  considered  -  for 
the  purposes  of  this  paper  -  to  be  the  baseline  against  which  the 
other  two  models  will  be  compared. 

Fig.  4  shows  the  polarization  curve  predicted  with  the  catalyst 
agglomerate  model  using  a  Nation™  film  thickness  of  S  =  80 
nm  and  demonstrates  that  the  model  is  capable  of  capturing  all 
the  relevant  features  of  a  typical  polarization  curve,  including 
the  sharp  drop  off  in  performance  at  high  current.  Three  load 
points  are  identified  on  the  figure,  corresponding  to  each  of  the 
‘classic’  regimes:  (A)  activation  dominated  losses  at  low  current, 
(B)  ohmic  dominated  losses  at  mid-range  current  and  (C)  mass 
transport  dominated  losses  at  high  current. 


-1  -0.5  o 
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interface.  Moving  down  the  channel,  the  location  of  the  max¬ 
imum  clearly  shifts  towards  the  centre  of  the  channel  region, 
indicating  the  increasing  importance  of  oxygen  mass  transport 
limitations  even  at  moderate  current  densities. 

The  sole  ‘free’  parameter  in  the  catalyst  agglomerate  model 
is  the  thickness  of  the  electrolyte  film  surrounding  the  catalyst 
agglomerates,  S.  This  value  was  estimated  to  be  approximately 
80  nm  in  the  two-dimensional  work  of  Sun  et  al.  [1]  in  part 
because  this  number  gave  reasonable  values  for  the  predicted 
mass  transport  limited  current  densities.  The  initial  three- 
dimensional  work  retained  this  thickness,  and  it  was  found  that 
the  predicted  current  density  dropped  sharply  owing  the  oxygen 
depletion  along  the  channel  which  is  not  accounted  for  in  two 
dimensional  models.  Accordingly,  this  parameter  was  varied, 
and  the  resulting  polarization  curves  are  depicted  in  Fig.  7.  As 
5  is  decreased,  the  limiting  current  increases  dramatically,  such 
that  a  value  closer  to  35  nm  is  needed  to  predict  overall  current 
densities  similar  to  the  two-dimensional  case  of  Sun  et  al. 


Current  Density  [mA  cm'2] 

Fig.  7.  Cathode  polarization  curve  for  varying  Nation™  thickness. 


Fig.  8.  Current  density  distribution  at  z/Dh  =  40  for  load  point  B. 


It  is  additionally  interesting  to  note  that  changing  the  elec¬ 
trolyte  thickness  also  changes  the  relative  importance  of  mass 
transport  losses  compared  to  electron  transport  (ohmic)  losses, 
as  evidenced  in  Fig.  8.  This  figure  compares  the  current  distri¬ 
bution  for  load  point  B  at  the  mid-point  between  the  inlet  and 
the  outlet  of  the  fuel  cell.  For  the  lowest  electrolyte  thickness, 
the  location  of  the  maximum  current  is  at  the  land/channel  inter¬ 
face,  with  a  significant  local  minimum  at  the  channel  centreline. 
However,  as  8  increases  the  current  distribution  begins  to  flatten 
almost  to  the  point  where  the  maximum  current  production  is 
located  under  the  channel  region.  These  results  emphasize  that 
there  is  clearly  a  significant  need  to  experimentally  characterize 
the  electrolyte  film  thickness  in  catalyst  layer;  particularly,  if 
accurate  predictions  of  current  gradients  and  overall  electrode 
performance  is  desired. 

4.2.  Comparison  of  catalyst  models 

The  vast  majority  of  (CFD-based)  fuel  cell  models  do  not 
yet  use  agglomerate  type  models,  and  as  such  represent  sim¬ 
plifications  with  respect  to  the  real  catalyst  morphology.  The 
agglomerate  model  is  itself  a  simplification  compared  to  the  true 
morphology  of  the  catalyst,  but  yet  less  so  compared  to  the  com¬ 
monly  employed  discrete- volume  or  thin-interface  approaches. 
Significant  mass  transport  losses  occur  due  to  dissolved  oxygen 
transport  in  the  electrolyte  phase  and  these  losses  are  accounted 
for  in  the  agglomerate  model.  In  the  absence  of  a  mechanism 
to  account  for  these  losses,  conventional  single-phase  catalyst 
models  are  unable  to  capture  the  sharp  drop  off  in  the  polar¬ 
ization  curve  observed  at  high  currents.  For  this  reason,  it  was 
decided  to  compare  the  results  of  the  agglomerate  model  to  the 
discrete  catalyst  and  thin-interface  catalyst  model  using  a  low 
film  thickness  of  10  nm.  This  tends  to  minimize  the  differences 
between  the  models,  and  as  such,  the  identified  differences  will 
be  more  significant  with  realistic  thicknesses. 

Fig.  9  shows  the  predicted  polarization  curve  using  each  of 
the  three  models.  In  the  case  of  the  discrete  and  thin-interface 
models;  the  parameter,  ascaie,  in  the  Butler- Volmer  equation  was 
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Fig.  10.  Current  density  distribution  near  the  inlet  for  the  three  catalyst  models 
(/  =  100  mA  cm-2). 


varied  such  that  the  curves  matched  at  very  low  currents  and  was 
then  held  constant  for  the  remainder  of  the  polarization  curve. 
This  method  of  altering  the  kinetics  is  common  when  comparing 
model  data  to  experimentally  measured  polarization  curves  for 
poorly  characterized  electrodes  (i.e.  Ref.  [9]);  but  is  unnecessary 
when  using  the  catalyst  agglomerate  model. 

It  is  not  surprising  that  the  thin-interface  model  predicts  the 
highest  potential  over  the  entire  polarization  curve,  as  it  does 
not  account  for  ohmic  losses  in  the  catalyst  layer  (ionic  or 
electronic),  nor  does  it  account  for  the  effect  of  the  reduced 
porosity  on  mass  transfer  in  the  catalyst  layer.  The  discrete- 
volume  catalyst  model  does  account  for  these  effects,  but  it 
does  not  account  for  oxygen  dissolution  in  the  electrolyte  phase, 
nor  does  it  account  for  the  effectiveness  of  catalyst  utiliza¬ 
tion  that  occurs  within  an  agglomerate  structure.  As  such,  the 
predicted  polarization  curve  for  the  discrete  model  closely  fol¬ 
lows  that  of  the  agglomerate  model  up  to  a  current  density  of 
approximately  200  mAcm-2.  It  then  exhibits  a  very  similar 
slope  to  the  thin-interface  model  until  a  current  density  of  about 
1400 mAcm-2  where  it  begins  to  capture  some  mass  transport 
losses. 

Since  both  of  the  differences  between  the  discrete-volume 
catalyst  model  and  the  catalyst  agglomerate  model  are  related 
to  mass  transport,  it  should  be  clear  from  this  figure  that  mass 
transport  effects  due  to  the  electrolyte  film  are  evident  in  the 
polarization  curve  at  current  densities  as  low  as  250  mA  cm-2 
(which  is  well  before  the  visible  inflection  point)  even  with  the 
use  of  a  minimal  electrolyte  film  thickness  of  10  nm. 

Figs.  10  and  11  compare  the  predicted  current  distributions 
near  the  inlet  and  near  the  outlet  of  the  fuel  cell,  respectively,  for 
each  of  the  three  catalyst  models  at  a  load  of  100  mA  cm-2.  The 
curves  are  all  normalized  by  iavg  =  100  mA  cm-2.  As  expected, 
the  discrete-volume  catalyst  model  and  the  catalyst  agglomerate 
model  predictions  are  very  similar.  The  thin-interface  model  is 
quite  different  however,  predicting  significantly  higher  currents 
near  the  inlet,  and  significantly  lower  currents  near  the  outlet. 
This  behaviour  stems  from  the  fact  that  the  kinetics  of  the  reac¬ 
tion  are  artificially  enhanced  through  the  scaling  of  the  effective 


catalyst  area  to  account  for  the  physical  morphology  of  the  cata¬ 
lyst  which  is  not  represented  at  all  in  the  thin-interface  approach. 
The  resulting  increased  consumption  of  oxygen  closer  to  the 
inlet  yields  reduced  kinetics  near  the  outlet  due  to  lower  oxygen 
concentrations.  The  net  effect  is  that  not  only  is  the  potential 
overpredicted,  but  current  gradients  are  grossly  over-predicted 
spatially  along  the  cell  using  the  thin-interface  model.  The  dif¬ 
ference  between  the  maximum  and  minimum  predicted  current 
is  0.48(avg  while  it  is  approximately  0.21  ;avg  for  the  other  two 
models. 

Fig.  12  compares  the  current  density  predictions  of  the  three 
models  near  the  fuel  cell  outlet  for  the  highest  load  case  of  zavg  = 
1500  mA  cm-2.  For  this  load  point,  the  thin-interface  model  and 
the  discrete- volume  model  predictions  are  similar,  with  the  thin- 
interface  prediction  consistently  lower  for  the  reasons  described 
previously.  The  catalyst  agglomerate  model  on  the  other  hand 
is  significantly  different  and  is  clearly  showing  evidence  of 
mass  transport  limitations,  consistent  with  the  physics  of  the 
approach. 
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Fig.  1 1 .  Current  density  distribution  near  the  outlet  for  the  three  catalyst  models 
(i  =  100  mAcm-2). 
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Fig.  12.  Current  density  distribution  near  the  outlet  for  the  three  catalyst  models 
(i'  =  1500  mA  cm“2). 

5.  Conclusions 

A  three-dimensional  implementation  of  a  catalyst  agglomer¬ 
ate  model  is  described,  and  applied  to  a  fuel  cell  with  straight 
flow  fields.  The  detailed  results  of  this  model  are  compared  to  the 
more  commonly  used  thin-interface  and  discrete-volume  cata¬ 
lyst  models.  It  was  found  that  the  catalyst  agglomerate  model 
is  capable  of  capturing  all  of  the  regions  of  the  polarization 
curve — including  the  sharp  drop  off  in  performance  at  high  cur¬ 
rents.  The  results  show  that  the  thin-interface  model  significantly 
over-predicts  the  current  gradients  across  the  fuel  cell  in  com¬ 
parison  to  the  two  other  approaches.  Parametrically,  the  catalyst 
agglomerate  model  also  shows  the  tendency  for  the  results  to  be 
quite  sensitive  to  the  single  remaining  ‘free  parameter’ ,  the  thick¬ 
ness  of  the  electrolyte  film  which  surrounds  the  agglomerates. 
This  elicits  the  need  for  strong  experimental  characterization 
of  the  film  thickness,  particularly  if  accurate  estimations  of  the 
mass  transport  related  losses  or  limiting  currents  are  desired. 

Finally,  even  though  relatively  high  stoichiometric  ratios 
were  used,  clear  evidence  of  mass  transport  ‘limitations’  are 
identified  at  very  low  currents,  suggesting  that  it  is  a  significant 
oversimplification  to  conceptualize  the  polarization  curve  using 
three  distinct  regions:  activation,  ohmic,  and  mass  transport.  This 
result  also  suggests  that  the  observed  experimental  mass  trans¬ 
port  limitations  of  electrodes  may  not  be  solely  attributable  to 
liquid  water,  but  that  the  electrode  structure  and  morphology 
may  also  play  as  significant  a  role. 
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